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Abstract 

The non-modal linear stability of miscible viscous fingering in a two dimensional homogeneous 
porous medium has been investigated. The linearized perturbed equations for Darcy’s law coupled 
with a convection-diffusion equation is discretized using finite difference method. The resultant 
initial value problem is solved by fourth order Runge-Kutta method, followed by a singular value 
decomposition of the propagator matrix. Particular attention is given to the transient behavior 
rather than the long-time behavior of eigenmodes predicted by the traditional modal analysis. The 
transient behaviors of the response to external excitations and the response to initial conditions 
are studied by examining the e—pseudospectra structures and the largest energy growth function. 
With the help of non-modal stability analysis we demonstrate that at early times the displacement 
flow is dominated by diffusion and the perturbations decay. At later times, when convection 
dominates diffusion, perturbations grow. Furthermore, we show that the dominant perturbation 
that experiences the maximum amplification within the linear regime lead to the transient growth. 
These two important features were previously unattainable in the existing linear stability methods 
for miscible viscous fingering. To explore the relevance of the optimal perturbation obtained from 
non-modal analysis, we performed direct numerical simulations using a highly accurate pseudo- 
spectral method. Further, a comparison of the present stability analysis with existing modal and 
initial value approach is also presented. It is shown that the non-modal stability results are in better 
agreement, than the other existing stability analyses, with those obtained from direct numerical 
simulations. 
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I. INTRODUCTION 


The Saffman-Taylor or viscous fingering (VF) instability, which occurs when a less viscous 
fluid is injected into a more viscous fluid, is a classical hydrodynamic instability phenomena 
[1, 21. Such instability occurs in many physical processes such as, enhanced oil recovery[I], 
spreading of surfactants coated thin liquid films [3], liquid chromatography [3], dispersion 
in aquifers [3j etc. VF has also been used as an effective mechanism to enhance the mixing 
of two fluids in a micro channel [5]. VF in miscible fluids were studied experimentally and 
theoretically by many researchers, such as Hill [6j, Slobod and Thomas [7j, Perkins et al. 
|8j, Tan and Homsy [9] are few to name. In recent times many substantial contributions 
ITOl [14] have been made in terms of experiments, developing new mathematical tools and 
numerical recipes for better understanding of this phenomenon. The main mathematical 
challenge arises in the linear stability analysis (LSA) of miscible VF due to the unsteady 
base state. Thus, to analyze its LSA one needs to solve a non-autonomous system of the 
linearized differential operators. For the time-periodic problems one can invoke the Floquet- 
theory [15] . But in an arbitrary time-dependent problem, Farrell and Ioannou [16] observed 
that the disturbances can no longer be set proportional to exp(crt), where cr is the growth 
rate. As a result, the linear stability theory of the non-autonomous system prohibits to use 
normal mode analyses. 

The following two approaches are common in practice for analyzing the linear stability of 
miscible VF to find the onset of fingering: (a) quasi-steady-state approximation (QSSA) and 
(6) initial value problem (IVP). In QSSA, one reduces to an autonomous system by freezing 
the base state at some specific time and apply modal analysis, while in IVP approach one 
can find the full solution of the non-autonomous linearized problem for some representative 
initial conditions. Tan and Homsy [9] reported that the QSSA is poorly agreeing with IVP 
and non-linear simulations. Although the IVP approach predicts the early time behavior 
better than QSSA, there are two major challenges in this method. The first one being 
how large the disturbances must grow before they become observable, and the other one 
is the choice of representative initial condition. The chosen random initial condition may 
not necessarily correspond to the one that gives optimum growth of the perturbation. Also, 
there is much disagreement on how to determine the growth of the perturbations from IVP 
analysis. Moreover, the random initial conditions, which are supposed to be localized within 
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the diffusive layer, perturb the system over the entire computational domain. 

Ben et al. |10j performed an LSA for VF using a spectral analysis method in a self¬ 
similar coordinate without invoking QSSA. Recently, Kim [12j analyzed VF in a miscible 
slice and compared the predicted growth rates from QSSA with those obtained using the 
spectral analysis. He found that the spectral analysis predicts the system to be initially 
unconditionally stable and becomes unstable at later times, while in contrary, QSSA [9] 
prediction reveals an initially unstable situation. As time progresses the predicted growth 
rates of the two methods get closer. A similar approach was also adopted by Pramanik and 
Mishra [H] to study the effect of the Korteweg stresses on miscible VF using a self-similar 
QSSA (SS-QSSA). 

However, the eigenmodes obtained from both SS-QSSA and QSSA are non-orthogonal, 
which indicates the possibility for transient growth of the disturbances mt It is a well- 
known fact that modal analysis does not address the phenomenon of transient growth in time 
na. In particular, it depends on the spectral properties of the underlying linear operator. 
For a non-normal operator (i.e., the operator that does not commute with its adjoint) 
an initial perturbation can grow in time by large factors before decaying, even if all the 
eigenmodes of the operator are damping. Earlier observations of transient growth due to the 
non-normality of the linearized operators were reported in the studies of instability in parallel 
shear flows | lR IT/] , many atmospheric and laboratory flows PH 120]. Most importantly, in 
all these hydrodynamic instability problems, the transient growth of infinitesimal small 
perturbations has been studied about a steady base state. 

In recent times, few studies pjj 12TH28] addressed the transient growth for flows with un¬ 
steady base state. Literature pertinent et al. [221123] used non-modal analysis (NMA) to 
determine the optimal perturbations with maximum amplification for density driven finger¬ 
ing by singular value decomposition (SVD) method. Later, Doumenc et al. [25] and Daniel 
et al. [28] used the ‘direct-adjoint looping’ method nznm for studying the transient growth 
of perturbations in Raylcigh-Bernard-Marogani convection and density driven fingering, re¬ 
spectively. Additionally, Daniel et al. [28] observed that the onset of convection in physical 
systems is due to suboptimal perturbations localized within the diffusive layer and hence 
proposed a modified optimization procedure. 

Despite the fact that transient growth analysis is known in various hydrodynamic insta¬ 
bility problems mentioned above, to the best of authors’ knowledge, the transient growth 
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of the perturbations is yet to be investigated for miscible VF. The non-autonomous linear 
equations and the non-orthogonal SS-QSSA eigenmodes, which may lead to the transient 
growth of perturbations, motivate us to pursue the non-modal analysis in such a system. 
In non-modal stability theory, stability is redefined in a broader sense as the response to 
general input variables, including initial conditions, impulsive and continuous external exci¬ 
tations [EJ [29]. In the present paper, the transient behavior and non-normality of the linear 
perturbed system have been investigated in terms of responses to continuous external exci¬ 
tations and to initial conditions. In order to know the response to external excitations, we 
computed e—pseudospectra and the study of transient growth of the optimal perturbations 
is performed through the propagator matrix approach. To quantify the maximum amplifi¬ 
cation of the perturbations we solve a matrix valued IVP, obtained from Darcy’s equation 
coupled with a convection-diffusion equation, by fourth-order explicit Runge-Kutta method 
to obtain the propagator matrix. Using SVD of the propagator matrix, the singular values 
and the right singular vectors are obtained. The obtained singular value and right singular 
vector provide the optimal amplification and the optimal initial conditions, respectively. At 
the early times, a substantial transient growth of the perturbations is observed to exist in 
miscible VF, which was not shown in the existing literature mmm of miscible VF. We 
also perform the direct numerical simulations (DNS) using a highly accurate Fourier pseudo- 
spectral method [30]. The results of the amplifications obtained by non-modal theory are in 
good agreement with those observed from the DNS. 


The paper is organized as follows. The governing equations along with appropriate bound¬ 
ary and initial conditions are presented in Sec. [TT] The linear perturbed equations and the 
numerical scheme are described in Section nm The rudiments of classical non-modal linear 


stability analysis is described in Sec. [IV] This part also describes the non-modal approach 
via singular value decomposition method. Section M discusses our numerical findings from 


NMA. Section [VT] contains the details about the DNS computations and non-linear growth 
rate. This section also discusses a comparison of previously studied linear stability methods 
with the results of NMA and DNS. 
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II. MATHEMATICAL FORMULATION 


We consider a uniform rectilinear displacement of a fluid of viscosity /12 by another 
fluid of viscosity n 1 bi a homogeneous porous medium, as shown in Fig. [lj The uniform 
displacement velocity is Ue x , where e x is the unit vector along the x— direction. The viscosity 
of the two fluids is determined in terms of a solute concentration, c, which is assumed to 
be c = 0 and c = c 2 in the displacing and displaced fluids, respectively. Fluids are assumed 
to be non-reactive, neutrally buoyant, and incompressible. The porous media have constant 
permeability, k, and an isotropic dispersion D. 



FIG. 1: Initially the interface is flat (dashed line) and then a wave like infinitesimal small 

perturbation is applied. 


A. Governing Equations 

The above-mentioned displacement flow in porous media can be described by the Dacry’s 
law for the fluid velocity coupled with a convection-diffusion equation that determines the 
evolution of the solute concentration. For the non-dimensional formulation of the problem 
we use U, D/U, D/U 2 , Hi, HiD/k as the characteristic velocity, length, time, viscosity, and 
pressure, respectively. The related non-dimensional equations in a reference frame moving 
with velocity Ue x are [9] 


<1 

IS 

0 

(1) 

Vp = -n(c)(u + e x ), 

(2) 

^+U -Vc = V 2 c, 

(3) 
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where u = (u, v) is the two-dimensional Darcy velocity. The only parameter that enters 
the dimensionless viscosity-concentration relation is the log-mobility ratio R, defined as, 
R = ln(// 2 ///i), and the viscosity is related to the concentration by an Arrhenius type 
relationship [9], i.e., //(c) = exp(ite). The coupled Eqs. ([lj - (J3| are associated with 
following initial and boundary conditions ED- 
Initial conditions: 


{ 0, x 0 

Wy. 

1, x > 0, 


Boundary conditions: 


I 

dc dv 

dy ’ dy 


= 0, |x| —> oo, streamwise direction, 

0, Wx, spanwise direction. 


(4) 

(5) 

( 6 ) 


B. Base State 

The base state of the flow is assumed to be pure diffusion of the concentration along the 
axial direction (i.e., no advection u = (0,0)). Hence, the time-dependent base-state solution 
of Eqs. ([Tj)- ([3]) on infinite domain with the no-flux boundary conditions for the concentration 

is i 


u 0 = (u 0 , Vo) = 0, 

(7) 

h-0 = A 4 o( c o) — ^o( x R), 

(8) 

p 0 (x,t) = - fjLo(s,t)ds, 

J —OO 

(9) 

co{x,t) = ^[l + erf(a:/2Vi)], 

(10) 


where erf(z) = f~e t2 dt is the error function. On transforming the current coordinates 
(x,y,t) to a new coordinate system (£,//, t), the base state concentration (Eq. (10)) becomes 


c o(0 = \ [1 + erf(f/2)], 


( 11 ) 


where ^ := x/y/t is the similarity variable transformation. Advantage of working with (£, y, t ) 
coordinate systems is that, one can pretend that as if the base state is time independent. 
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III. LINEAR PERTURBATION EQUATIONS 


In this section, we discuss the method of numerical solutions of the linearized equations. 
The linear perturbation equations are solved as an IVP in both (x, t) and (£, t) coordinate 
systems, and the obtained numerical results are compared. The response of the system to 
small disturbances is obtained by linearizing Eqs. 0-0 according to 

, V , t) = u'(£, y, t ), c(f, y , t) = c 0 (£) + c'(£, y , t), (12) 


where u'(£,y,t), d(£,y,t ) represent infinitesimal disturbances to the velocity and concen¬ 
tration, respectively. We seek disturbances of the form ( u',d)(£,y,t ) = («', <:')(£, t)e lky , 
describing a function which evolves in time and in the streamwise direction while exhibiting 
sinusoidal character in the transverse direction. Using viscosity-concentration relationship 
//(c) = e Rc , the disturbances u' and d are determined from the following set of equations 
representing a disturbance of transverse wave number k 


V 


R 

20F 


exp 




d_ 


d £8 _ 

t - X> 

dt 2 8^ 


4 

c'(eu) = 


u' (£,t) = k 2 Rtc'(£,t ), 


Vi 


2VH 


exp 


-e 2 




(13) 

(14) 


where V := d 2 /d£ 2 — tk 2 . The associated boundary conditions are ( c',u ') —> (0,0), as 
|£| — y oo. Note that u' and d are not periodic in the ‘£’ variable due to the far held 
conditions i.e., ( d,u') — >■ (0,0), |^| — > oo. So the functions v! and d cannot be represented 
as a superposition of sinusoidal perturbations in the ‘£’ variable. The analytic solution 


for the coupled system of partial differential Eqs. (13) and (14), for a diffusive base state 
is unattainable. Hence, a numerical solution method needs to be used for solving these 
coupled equations. In this paper we use finite difference method and solved the linear 
stability problem in a finite computational domain. 


A. Numerical Method 


For a given k and R, we discretize Eqs. (13) and (14) over a finite computational domain 


[—L,L\, L e M + , with a uniform mesh size h and n + 2 number of grid points. All the 
derivatives are approximated by central difference formulae. The discretized version of Eqs. 
(13) and (14) with the boundary conditions u'(t) = 0 = d(t) and u'(t) 


i=0 


i=0 


i=n+l 


7 




















FIG. 2: Dependence on the size of the computational domain on growth rate obtained from 
NMA, for R = 3, k — 0.2. Beyond the domain length L = 90, there are almost no change. 


0, c\t) 


= 1 can be cast into a non-autonomous system 

i=n -\-1 


d c'{t) 
df 


A(t)c\ 


(15) 


where c'(f) = c'(£i,£) and A(t) = M 3 + is the stability matrix of order n. Here 

each matrix Mj , j = 1,2, 3,4 is either a diagonal matrix or a tridiagonal matrix and are 
given by [Q], 


p ± ( c o U + !) - c 0 (j - 1)), i-jTl 


Mi (i,j) = < 


_ U2 f 

h 2 F 

0, 


* — J 

otherwise, 


and 


M 2 (i,j) 


{ ik 2 Rt , i = j 
0, otherwise, 


M 3 (i,j) = < 


th 2 =*= 4ht ’ ^ — J =F 1 

0, otherwise, 


M^iJ) = 


co(j+l)-co(j-l) 

2hVt 


1 = 2 


0, 


otherwise. 




















For the computational purpose the unstable interface is set at £ = 0. The initial value 


problem, Eq. (15), is solved using Runge-Kutta fourth-order explicit method for time in¬ 
tegration with a relative error of order (9(10“ 10 ). The matlab built-in functions ode45 and 
eigs are used for solving the IVP and Ending the eigenvalues, respectively. 

The computational domain is taken sufficiently large so that the numerical results remain 
unaffected by the truncation of infinite physical domain into finite domain in numerical 
algorithms. To ensure that our numerical quantities, viz. singular values, eigenvalues, etc., 
that characterize the instability, are independent of spatial domain a series of numerical 
calculations are performed for different L ranging from 50 to 120. The parameter values 
are chosen to be R = 3, k = 0.2 [9]. The obtained numerical results are compared in Fig|2j 
and it is observed that for L > 90 there is no change in the growth rates calculated from 


NMA which is explained in Sec. IV B Thus for optimal result, the computational domain 
is chosen [—100,100] in all our calculations. 



-10 -5 0 5 10 

£ 


(b) 



FIG. 3: Disturbance concentration prohle for R = 3, k = 0 at different times, obtained 
from the IVP in (a) (£, t) coordinate system, (6) ( x,t ) coordinate system. 


Numerical convergence of finite difference method has been tested with different number 
of grid points, and the optimal computational time and accuracy were obtained with the 
spatial step size 0.2. As a verification of our numerical code we reproduce the results of 
Rapaka et al. [22]. Due to the possible singularity at t — 0, all the numerical integrations 
are performed starting from the initial time t = 10~ 3 to various final time. Although we 
restrict our study of the optimal perturbations and maximum amplification for R = 3, the 
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analogous transient growth dynamics of perturbations can be observed for other positive 


values of R. The next section (Sec. Ill B) describes the importance of NMA for the present 


problem by using the non-autonomous system (Eq. (15)) and explains the result of IVP 
calculations in (x, t ) and (£, t) co-ordinate systems. 


B. Initial Value Calculation 


In this section, we illustrate the advantages of studying stability analysis in (£, t) domain. 


For this, we solve the non-autonomous system Eq. (15) with respect to the following random 
initial condition 


c'(uj, t 0 ) = 5 * rand (a;), where lu = £ or x, (16) 

where rand(cu) represents a random number generator between —1 and 1, and an infinitesimal 
parameter 6 determines the amplitude of the perturbation, which is chosen to be (9(10~ 3 ). 
The objective of performing IVP in both the coordinate systems is to describe the structure 
of perturbed concentration in each coordinate system. The initial time to = 10~ 3 has been 
chosen for comparing the structure of perturbed concentration obtained from (x, t ) and 
(£, t) coordinates. Fig. J3|a) depicts that the concentration eigen-function determined by 
solving the IVP in (£, t) domain converges rapidly to the dominant eigenmode exp(—£ 2 /4) 
[TO] . On the other hand, Fig. |3](b) illustrates that the concentration eigen-function obtained 
from the IVP calculation in (x,t) coordinate system takes longer time to converge to the 
dominant mode, in comparison to calculations in (£, t) coordinate system. Both the results 
shown in Fig. [3] are obtained by averaging over 10 random initial conditions. The localized 
eigen-functions in (£, t) coordinates converge rapidly to the exact solution, thereby giving an 
accurate disturbance growth rate at small times. With this observation, rest of the numerical 
calculations are performed in the self-similar (£, t) domain. We perform IVP calculations in 
(£, t) domain with dominant eigenmode exp(—£ 2 /4) as the initial condition and the obtained 
results are discussed in section IVI A II 


IV. TRANSIENT BEHAVIORS AND NON-MODAL INSTABILITIES 

In this section, the stability analysis is carried in a broad sense as the responses to 
external excitations and to the initial conditions. We describe both of these through rigorous 
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mathematical analysis. Suppose the system is subjected to a perturbation of the form 


«'(£,*) = ii'K)exp(cr(«o)(), c'K,«) = c'(f) exp(<r(t 0 )t),, 


( 17 ) 


where a is the complex frequency, and to is the diffusive time at which the base state is 


frozen. Thus Eq. (15) is converted into an eigenvalue problem A(t 0 )c. = a(t 0 )c'. To analyse 
the system behavior to the external excitation the study of e—pseudospectra[I7, 23 SZJ is a 


very helpful tool. For given e > 0, the e—pseudospectra is defined as 


K(A) = {* 6 c :|| RM) ll>« -1 }. 


(18) 


where R Z (A) = (zl — A)~ is known as the resolvent, I is the identity operator, and || • || 
is the Euclidean 2—norm in C n . 

In pursuit of optimal amplihcation of initial conditions, we follow the classical propagator 
matrix approach. The maximal amplihcation is obtained as the singular value decomposi¬ 


tion of this propagator matrix and the detail analysis is given in Sec. IV B This section also 
describes the choice of measure in terms of Euclidean 2—norm. Often for IVP in hydrody¬ 
namics, one may be interested in the dynamics of the energy growth rate, for t -C 1. In this 
context, the concept of numerical range can be used to link the stability matrix to the early 
time energy growth H3E2]. The numerical range is defined as 


W(A) = {xMx : x G 


x 11= !}. 


(19) 


where the superscript * denotes the transpose. The growth or decay of the initial energy 
can be studied from the numerical abscissa 


rj(A) — sup = sup{A(A + A*)/2}, (20) 

zeW(A) 

where A(-) represents the spectrum of the matrix and denotes the real part. The numer¬ 
ical abscissa r](A) measures the maximum possible instantaneous growth rate corresponding 
to any initial condition as t —> 0 [32] . 


A. Pseudospectra 

The behavior of the eigenvalues of the linearized stability matrix A(t) determine the in¬ 
stability of the system. But, for the highly non-normal matrices, e.g., in bounded shear flows 
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(a) (b) 



(c) (d) 




FIG. 4: Spectral abscissa (circles) and numerical abscissa (squares) for k = 0.2, R — 3, at 
different time (a) t = 10~ 8 , ( b ) t = 10 -5 , (c) t = 1, (d) t = 10 . Also shown the eigenvalue 
with largest real part (black dot) and imaginary axis as bold continuous vertical line. 3ft 
and A denote the real and imaginary axis in the complex plane. 


PH EH, the classical normal mode analysis, which assumes an exponential time dependence 
of A(t), fails to predict the instability appropriately. The temporal eigenmodes of the time 
dependent matrix A(t) are not defined [TB] and to determine the asymptotic behavior (by 


the Lyapunov exponents) of Eq. (15) is beyond the scope of this paper. To quantify the non¬ 
normality of A(t), we freeze A(t) at different times and consider two quantities namely the 
spectral abscissa and the numerical abscissa. The spectral abscissa of the stability matrix 
A is defined as 


a(A) : = max{3ft(A(A))} (21) 

It can be verified that for a normal matrix a(A) = The main objective of modal 

analysis is to study the spectral abscissa and the corresponding eigenmodes. In Fig. [4| the 
spectral and numerical abscissa for the stability matrix A are plotted at different time with 
for given wave number k = 0.2 and R — 3. It is seen from Figs. S») and [4|b) that at 
early times, the spectral abscissa and the numerical abscissa differ from each other in the 
order of 0(10') and 0(1O 4 ). This reveals that at early time, the stability matrix A is highly 
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(a) 


(b) 




FIG. 5: Pseudospectra of A(t) for (a) R = 2, k = 0.15, £ = 5, (b) R — 3, k — 0.25, £ = 3. 
Black square: numerical abscissa; black dots: eigenvalues; dashed line: boundary of the 
numerical range; solid lines: contours from innermost to outermost representing levels from 
e = 10 -2 to 10 -0 ' 5 with increment lO 0,5 . Here 3ft and S represent the real and imaginary 

axis in the complex plane. 


non-normal, but at later times, A(t) is close to a normal matrix (see Figs. 0c) and[f](d)). 


For the unsteady base state, Pramanik et.al studied the response to external excitations 
by analysing the pseudospectrum [33]. Following these authors, we plot e—pseudospectra 
using Eigtool [33]. Fig. [5] illustrates the pseudospectra for R = 2, k = 0.15, £ = 5 and 
R = 3, k = 0.25, £ = 3. The contour lines indicate isolines of constant resolvent norm. We 
observe that the numerical range partially reaches into the unstable half-plane and it strictly 
contains the spectrum. This means that there exists positive energy growth rates, despite 
the fact that the spectrum is contained inside the stable half-plane. This also shows that 
the eigenvalues are inherently time-asymptotic tools when dealing with stability analysis of 
non-normal matrices. We thus conclude that initial energy growth is possible up to a growth 
rate given by the maximum protrusion of the numerical range into the unstable half-plane 
and the least stable eigenvalue governs the long time behavior. This early time structure of 
non-normality of A(t) inspires us to investigate the transient growth of perturbations and 


the onset of instability, which is presented in Sec. IV B 
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B. Transient energy growth 


In this section, we present the mathematical analysis for calculating the transient energy 


growth and the onset of instability. Let us assume a solution of Eq. (15) as 


c'(f) := <S(«o;i) c o. 


( 22 ) 


where c'{t Q ) = c' 0 is an arbitrary initial condition. Here Q(t 0 ]t) is called the propagator 
matrix, because it propagates the information forward from the initial time t 0 to time t. On 


substituting Eq. (22) into Eq. (15) we obtain a matrix differential equation, 


—#(«„; i) = A(t)$(t 0 ;t), 


(23) 


with the initial condition <E>(£ 0 ; to) = /, where / is the identity matrix of the same order of 
the matrix A(t). Thus, instead of solving a vector differential equation with random initial 


condition, now it needs to solve a matrix differential equation, Eq. (23), with deterministic 
initial condition, i.e., <E>(£o;to) = I- 

The stability analysis is performed based on the amplification magnitude of the perturba¬ 
tions over a prescribed finite time interval HU- In the stability analysis of a physical system 
with time dependent base state, the growth or decay of disturbances is only meaningful in 
reference to the evolved base state. In 1961, Shen m observed in his study of unsteady 
parallel shear flow that if both the disturbance and the base state are decaying, but the 
latter one with faster rate, the disturbance, relative to the basic state, would appear to be 
amplified at a later instant. Conversely, if a disturbance grows in time, but the base state 
grows faster than that, then the disturbance will appear to decay in time. So it is necessary 
that the growth or decay of a disturbance should be measured only by a comparison with the 
basic flow. With these observations Shen EH introduced “momentary stability” and defined 
an appropriate measure for stability analysis. Similar measures have also been adopted by 
Matar and Troian [3] in their studies of spreading of surfactants on a thin film, and Bestehorn 
and Firoozabadi (26] in their study of the dissolution of C0 2 in saline aquifers. 

Following the work of Shen EH, we would like to introduce the quantities such as amplifi¬ 
cation magnitude and transient growth rate of the disturbances that will be used to measure 
the time dependent growth rate of the perturbations. The perturbation magnitude for the 
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flow variable / can be defined as 


M,(t) := 


E,(t) 

Efjt)’ 


(24) 


where f 0 = c 0 ,or u 0 , and / = c', or u' and E f (t) - ||/(()lll = P(x,t)dx, E h (t) - 
|/ij!di ||) — J_ /q (x, f)dx. The sensitivity of the infinitesimal disturbances introduced at 
time to is measured from the ratio of the perturbation magnitude Mf(t ) at time t to its 
initial value Mf{to)- Consider the normalized amplification T/(t) defined by 


M,(t) 

E,(t) 1 

/ 

E f0 ( t ) 

G f (t ) 

Mf(t 0 ) 


/ 

.Ef 0 (to). 

G f0 (t) 


*f(t) ■ = 

The time-dependent growth rate is defined as |3], [2HJ 

a(t) - 1 d'I' / (t) _ 1 d Gf 1 d G fn 


(25) 


= a/(f) - a h (t). 


(26) 


T/(t) dt Gf dt Gf 0 dt 
As Co in (£,f) coordinate is independent of time and uq = 0, this implies cr/ 0 (f) = 0, and Eq. 


(26) reduces to 


a(t) = 


1 d G f 


(27) 


Gf d t 

To measure the transient growth, we consider the concentration perturbation magnitude, 
i.e. / = d at time t which is normalized to the concentration perturbation magnitude at to- 
As the perturbation equations are linear, without loss of generality it can be assumed that 


|c'(to)|| = 1 and using Eq. (22), Eq. (25) reduces to 


T c /(t) = E c ,{t)/E c ,{to) = ||c'(a:,t)|| 2 
= (§(to\ty<$>(to\t)c! ,c'), 


(28) 


where (-,-) is the standard L 2 inner product. The objective of NMA is to find the largest 
amplification, G c /(t) = G(t), say, which the disturbances can achieve over all possible initial 
conditions, i.e., 


G{t) = G(t,k,R ) := maxT c /(t) = max ||$(t 0 ; t)cg ||2 

c 0 c 0 

= ll $ (^o;t)|| =supsj(t), (29) 

3 

where s/s are the singular values of <h(t 0 ;t), that is the eigenvalues of $(t 0 ; t)*$(t 0 ; t) and 
can be found by SVD of <h(to; t). Thus, the quantities that determine the optimal amplifica¬ 
tion and the finite time behavior of energy growth are eigenvalues of <I>(£o; t)*$(to; t). It must 


15 


















be noted that we do not calculate the E c /(t) explicitly, instead, we are using the classical 
definition of amplification magnitude of disturbances to find the optimum amplification and 
the optimum perturbations by SVD. Further, there exists other possible ways to describe 
such transient growth instability [241128 ]. For example, in the study of densdity-stratifed flow 
Daniel et al. |28j modified the definition of Ef(t) (given in Eq. ([24]) ) to accommodate both 
concentration and velocity disturbances i.e., E e (t) := [c /2 (a:, t) + u' 2 (x, i)] da: (see Eq. 

(3.1) of Daniel et al. [28]). It was concluded that for such flow problem, the amplification 
measure E c /(t) was sufficient for their stability analysis. In an another study of Rayleigh- 
Bernard-Maragoni convection, Doumenc et al. [25] observed that the velocity perturbations 
are slaved to the temperature perturbation when the time derivative of velocity was dropped 
out from their model equations. Hence they discussed the stability criteria based on the op¬ 
timal amplification associated to the temperature disturbances only. Similar observation is 
also noted by Ward et al. [35] that the stream function is slaved to the concentration field, 
in the study of convection of reactive solute in porous media. Thus, following the above 
observations and the face that the linearized equations (Eqs. (13) and (fl4|)) do not contain 
the explicit time derivative of velocity, the concentration disturbances are pertinent in the 
present problem. Hence, the optimal amplification associated to the concentration distur¬ 
bances is used for the non-modal stability analysis. Below we define the growth rate and 
the critical time for the quantitative investigation of the transient growth of perturbations. 


Definition IV. 1 (Stability criterion and Growth rate). The instantaneous growth rate a 
at time t of a mode with wave number k can be defined as: 

=a{R ’ k ’ t):= W)^r" (30) 

A given perturbation is stable if a(t) < 0 and unstable for cr(t) > 0. With this definition 
of growth rate the critical time or the onset of instability for the given perturbation can be 
found. 


Definition IV.2 (Critical time). For a given perturbation the critical time is denoted by 
t c (k,R ) and can be defined as the first instance of time when 


= 0 . 


t=t c 


As noted by Doumenc et a/.[25j that the optimum amplification (see Eq. (29)) has to be 
determined with subject to the following constraints: (a) the disturbance energy at initial 
time t 0 is equal to unity, (b) the disturbance must satisfy the linear governing equations 
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(a) (b) 




FIG. 6: (a) Optimum amplification in (x, t ) (dashed line) and (£, t) (continuous line) 
coordinate systems with k = 0.2, and R = 3. In (£,t) coordinate system, the optimal 
amplification log 10 (G(t)) is a non-monotonic function of time, which signifies the 
dominance of diffusion at early times. Also shown the onset of instability (black dot) 
determined by non-modal analysis in (£, t) domain, (b) Optimal amplification in (£, t) 
coordinate system, for R = 3 and fc = 0.01, 0.1, 0.3. 


along with the boundary conditions in the time interval [t 0 ,t]. With this observation i.e., 
G(to) = || < h(t 0 ;O)|| = 1, a relationship between the optimum amplification, G(t), and the 


growth rate, cr(i), can be obtained from Eq. (30) as [27] 


G(t) = exp 


cr(s)ds 


'to 


(31) 


V. RESULTS AND DISCUSSION 

First we explain the findings of NMA, which includes the optimal amplification of the 
disturbances and the spatial structure of it. We illustrate that the optimal perturbations 
are localized about the base state, which is contrary to the random initial perturbations 
that perturbs the system across the entire domain. Then, we analyse the onset of instability 
obtained from DNS, NMA, and other existing LSA. 
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A. Amplification 


In this section, the optimal amplification G(t) of a given perturbation measured from 
(x,t) and (£, t) coordinate systems by NMA is presented. For the base state cq(£, t) (see 


Eq. (11)), given a value of log-mobility ratio R, time t and wave number k, the non-modal 
analysis determines the optimum amplification over all possible perturbations. We choose 
a set of wave numbers k = 0.01,0.1, 0.2, 0.3 and plotted the optimal amplification G(t) as a 
function of time in Fig. |6j 

For k = 0.2 and R = 3, Fig. [6]( a) illustrates the optimal amplification log 10 (G(t)) obtained 
both in (x, t) and (£, t) coordinate systems. It is observed that the optimal amplification 
is a non-monotonic function of time in (£, t) coordinate system. The reason for this non¬ 
monotonicity is that at early times, the effect of diffusion is prominent which causes all the 
perturbations to decay. But at later time when this diffusion is weaken and convection is 
dominated, we observe an increase in amplification. The time when G(t) starts increasing 
is considered as the onset of instability, which is embedded in our definition of instability 
criterion definition IV.2, Further, from Fig. |6](a) it is observed that after the onset time G(t) 
reaches a transient peak, which at later time must follow an exponential decay. This physical 
characterization of the system is appropriately captured by NMA in (£, t) coordinate system. 
This is a possible reason for linear stability analysis in the self-similar coordinate system 
(£, t) to give a reasonable explanation to the actual physical situation. Since the definition of 


amplification G(t) (see Eq.(29)) is obtained by taking all admissible perturbations, it include 
the mode that provides maximum growth in SS-QSSA. Thus, the curve G{t) vs t can be 
thought of as an envelope over optimal initial conditions. In Fig. |6](b), it is shown that for 
R = 3, k = 0.1 and k = 0.3, G(t) is a non-monotonic function of time, except for k = 0.01, 
for which G(t) decays monotonically. Thus, there exists a range of wave numbers for which 
transient growth is prominent. This transient growth of perturbations is consistent with our 


observations in Sec. IV A obtained in the study of numerical and spectral abscissa. In that 
section we have shown that at early time there is a substantial difference between a(A) and 

v(A)- 


We observe that the range of wave number between 0.1 and 0.3 dominate most, and the 
most unstable wave number predicted from the modal analysis [9J lies within this interval. 
Thus, it can be claimed that the non-modal theory not only determines the transient growth, 
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FIG. 7: Contour plot of optimum initial perturbation at time t = 0.1, k = 0.2, and R = 3 
(dotted line correspond to negative contour lines and continuous line corresponds to 
positive contour lines). Inset shows the base state concentration cq. 


(a) 



£ 


FIG. 8: Optimal normalized perturbation V 0 


(b) 



£ 


(a) for k = 0.2 and R = 3 at different times, 


with the eigenfunction calculated from SS-QSSA at t — 3.5, ( b ) for R = 3 and 


k = 0.15,0.175,0.2,0.225 at their respective critical times calculated from definition 


IV.2 


but also is a generalization of modal and IVP analyses. It is also observed when k —y 0 (zero), 

Eqs. (13) and (14) with boundary conditions Eqs. (J5]) and @ simplify to 

d £ d d 2 ' 
dt 2 <9£ 2 


c'(£,t) = 0. 


(32) 


From Fig. [6](b) it is shown that unlike density fingering[28], the amplification is a non¬ 
constant function of time for any finite wave number. This is in good agreement with 
the previously studied linear stability analysis of VF in the self-similar coordinate systems 


mmm- 
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B. Optimal Perturbations 


One of the most important aspects of the non-modal analysis is to investigate the distur¬ 
bances that produce maximum response during the fingering process. It is not clear from the 
existing linear stability analyses of miscible VF what will be the dominant perturbations at 
early time on the order of t ~ t c . This section provides physical explanation of the transient 
period by examining the optimal perturbations and eigenmodes obtained from SS-QSSA. 
The maximal amplification G(t ) = ||<3>(to,t)ll can be found by computing its singular value 
decomposition or Schmidt decomposition [5E] , 


<F(f 0 ,t) = U [M^ [ t ot ]V’Wb ( 33 ) 

where V* represents the adjoint of the matrix V. The largest singular value is the optimal 
amplification for the given perturbation i.e., the first entry of the diagonal matrix Xl[t 0 tp 
and the corresponding column of V (first right singular vector) is the optimal perturbation 
(initial condition), denoted by V opt . This optimal initial condition evolves into the first 
column of U (the first left singular vector), denoted by U opt . This is also evident from 
the singular value decomposition of <f>(f 0 ,Q, i.e., <h(fo, t)v\ = Sgiib, where s^i is the first 
entry of the diagonal matrix )T )u t ,. Thus, the leading vectors of V are the useful inputs 
for selecting the correct representative for optimal initial conditions. Fig. [7] shows the 
contour of the initial perturbation of concentration that is most amplified at t = 0.125. 
As expected, the perturbation is localized within the diffusive layer. In Fig. |8](a) we plot 
the normalized optimal perturbation V opt for fixed k = 0.2 at time t = 0.125, 0.25,1, 2, and 
4.125. It is shown that these right singular vectors are situated around the base state and 
gradually resolves to the dominant eigenmode calculated from modal analysis discussed by 
Ben et al. [10]. In Fig. [8](b) the normalized V op t, for various wave numbers are shown at 
the respective critical time (see definition IV.2). Two important observations are noted 
from the structure of the optimal perturbations determined from non-modal theory: (i) 
the optimal perturbations at their respective critical times are identical to the eigenmode 
predicted by SS-QSSA, irrespective of any wave number chosen, and, (ii) at early times the 
SS-QSSA eigenmode and the optimal perturbations differ substantially, which is not shown 
for brevity. This reflects that although SS-QSSA captures the early times behavior better 
than QSSA and predicts the onset of instability, it fails to capture the transient growth of 
perturbations. Hence, a physically relevant transient period for VF can be defined as [t 0 ,t c ], 
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where to is the time when the time integration starts and t c is the critical time obtained 
from NMA. 


VI. NONLINEAR SIMULATIONS 


In this section, we discuss the direct numerical simulations for solving the coupled non¬ 
linear equations that are compared with the transient growth obtained from NMA. We solve 
the nonlinear problem using a highly accurate pseudo-spectral method based on stream func¬ 
tion formulation of Eqs. 0-0, proposed by Tan and Hornsy [30]. This method has been 
employed successfully to obtain highly accurate results for various viscous fingering problems 
[4], Writing the unknown variables as, y, t ) — ipo(x, t ) = ip'(x, y, t ), c(x, y, t ) — co(x, t ) = 

c'(o;, 2 /, t), the stream function form of Eqs. (jl|-([3]) in terms of the perturbation quantities 
can be represented by, 


vy = -r (W • v (cq + d) + , 


dc' d'd' (dc 0 dd\ 


dt dy 


dx dx ) 


W 9d _ V 2 c , 

dx dy 


(34) 

(35) 


Here the base-state flow is the same as discussed in Sec. |IIB[ i.e. yy = 0, c 0 is given by Eq. 
(J7]) and the primes correspond to their perturbation quantities (not necessarily infinitesimal). 
The non-dimensional width of the computational domain becomes Pe = UH/D , the Peclet 
number and the corresponding length of the domain is A ■ Pe. Here A = L/H is the 
aspect ratio and L , H being the dimensional length and width of the computational domain, 
respectively. Periodic boundary conditions are applied in the transverse direction. In the 
longitudinal direction, we work with the periodic extension of a displacement front. In this 
paper, the computational domain is chosen in such a way that Pe = 512 and A = 8 with 
1024 x 128 grid points discretizing the domain. The time integration is performed by taking 
time stepping 10~ 3 . 

The simulations are performed with to = 0.01 and initial concentration perturbation of 
the form, 


d(x, y, 0) = e rand(-) sin iky). 


(36) 


Here, k is the wave number, rand(-) represents a random number between 0 and 1 at the 
diffusive interface, e is the amplitude of the perturbation, which is taken to be 10~ 2 . As¬ 
suming that the perturbations grow exponentially, the growth rate of perturbations, (Tons, 
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can be computed at every time from, 


_ 1 d[ln(£ , DAfs(^))] 

aDNS = 2 df 


(37) 


Here, E DNS (t ) = fp[c'( x,y,t)] 2 dxdy represents an amplification measure of the per¬ 

turbed concentration. The DNS results are averaged over ten different random realizations 


of the initial concentration perturbations, Eq. (36). For our analysis DNS are performed for 
seven different wave numbers and the results obtained are discussed in Sec. IVI Al 


A. Comparison of non-modal theory with SS-QSSA, IVP and DNS 


In the previous sections (sec. V A and VB), it is shown that the structure of the optimal 
perturbation describes about the regime of transient growth. Also from the discussion of 


numerical and spectral-abscissa (see Sec. IVA), we noted that at early times there exists 
significant difference between a(A(t)) and r] (A(t)), which is of order O(10‘) (see Fig. [4j. To 
explore the relevance of this optimal perturbation to the physical system, we compare the 
results of non-modal analysis with those obtained from SS-QSSA, IVP and DNS. 



FIG. 9: Comparison of dispersion curves obtained from SS-QSSA (line with asterisk), IVP 
(line with square), NMA (continuous line), and DNS (circles) for R = 3 at time t = 3.5. 


1. Dispersion Curves and Growth rates 

Fig. |] shows the dispersion curves obtained from non-modal analysis, SS-QSSA, IVP 
in (£,f) domain, and DNS at t = 3.5. The relevance of plotting the dispersion curves at 
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t = 3.5 is that it is close to the onset of fingering obtained from non-modal analysis (see 
Fig. 11). It is to be noted that in Fig. [9] (also in Fig. [XT]) the time for SS-QSSA calculation 
must be interpreted as diffusion time, where the base is frozen. It is evident from Fig. 
[9] that non-modal analysis captures the onset of the instability more accurately than any 
other linear stability methods. Although both IVP and SS-QSSA show the flow is stable, 
the non-modal analysis is meticulously agreeing to DNS calculation. Moreover, it is shown 
that the threshold wave number (the least wave number at which a = 0) and the cut-off 
wave number k c (wave number for which a changes from positive to negative) obtained from 
non-modal analysis are closer to DNS than those obtained from SS-QSSA and IVP. Thus, 
there exists notable transient effect of linearized non-normal matrices in VF and the present 
nonlinear simulations agree with the non-modal linear stability results at early times. 



FIG. 10: Comparison of the amplifications from non-modal theory (continuous line) and 

DNS (dashed line) for R = 3 and k = 0.25. 


With these agreement of growth rate in NMA and DNS, we plot the amplification obtained 
from DNS and non-modal theory in FigjlOjfor the wave number k = 0.25, which amplifies 


most at early time. The amplification for DNS calculations are obtained using Eq. (31). 
It is shown that non-modal theory has an excellent agreement with DNS calculations up 
to f ~ 15. This confirms the physical relevance of the transient growth calculation from 


non-modal theory and the instance, when the two curves in Fig. 10 deviate from each other 
can be marked as the threshold of nonlinear convection. 
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0.5 


(a) 


(b) 



FIG. 11: (a) Optimum growth a max t over all possible wave numbers k, for R = 3, (6) 
Neutral stability curve comparison for R = 3: Continuous lines for SS-QSSA and dashed 

lines for NMA. 

2. Optimal growth and dominant wave numbers 

In an experiment or in any physical process a perturbation consists of combination of 
different wave numbers. Therefore, it is important to calculate the onset of instability by 
considering all possible wave numbers. This will characterize what will be the optimum 
growth at a given time. Fig. [fjja) depicts that the onset time obtained from modal, non- 
modal, and DNS are 4.6, 3.2, and 2.8, respectively. One of the important observations from 
optimal growth calculation is that the effect of the diffusion at early time is well described by 
non-modal and DNS, whereas SS-QSSA fails to observe this physical phenomenon effectively. 
The reason that NMA explains the physical phenomenon appropriately is attributed to 
the fact that at early times the linearized matrix A(t) is highly non-normal, which is not 
reported in earlier linear stability analyses. Since at very early stage of the evolution of the 
perturbations, the growth rate in NMA and SS-QSSA are all negative, it does not affect the 
instability. 

Thus, at early stage (i<Cl) in the linear regime optimal growth is damped before it starts 
to grow. The instance when both the convection and the diffusion terms are balanced, i.e. 
when a = 0, is termed as the neutral stability condition. It is often required and useful to 
identify the unstable and stable regimes precisely. Fig. [TT^b) represents the neutral stability 
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curve determined by modal and non-modal theories. Again, the qualitative features remain 
the same, but there is a quantitative difference at early times and hence the stable regime 
identified by modal analysis becomes unstable for NMA. The lowest point on the curves 
corresponds to the critical time and critical wave number. This clearly depicts that although 
the critical wave numbers obtained from both modal and non-modal analyses are the same, 
the critical time is smaller in non-modal theory than modal analysis. 


VII. CONCLUSIONS 

A novel non-modal linear stability analysis is presented to study the transient growth of 
perturbation in miscible VF. Recent works|2l L22, [28] suggested that in the unsteady base 
flow problem, the transient response of the perturbations can be significant. The existing 
literature does not discuss the transient effect of time dependent base state in the case of 
viscosity unstable problems. Our approach of non-modal analysis of VF is based on finding 
the propagator matrix and singular value decomposition method [22, [23]. This method 
predicts the dominant perturbations that experience the maximum amplification within 
the linear regime. The present linear stability analysis based on sounder mathematical 
basis have two most important features: (i) it determines the optimal amplification of the 
perturbation that is imperceptible in both SS-QSSA and IVP analyses, and (ii) it identifies 
the physical mechanism of the VF instability, which is in very good agreement with DNS 
results. Additionally, the structure of the initial perturbations which lead to the optimal 
amplification is also illustrated, along with all the qualitative information of flow instability 
that SS-QSSA and IVP together provides. Unlike the random perturbations which disturb 
the system everywhere, the initial perturbations obtained from NMA are localized within the 
diffusive zone, which is in agreement with the spectral analysis of Ben et al. [10]. Although 
matrix differential method is used in the present study, the adjoint-looping method [20, [25], 
[25] can be an useful alternate way to study the transient growth. Further, this classical non- 
modal analysis can be helpful to find the suitable time and length scales for the instability, 
the structure of most amplified initial condition in many important hydrodynamic instability 
problems, which deals with the unsteady base state. For example, the effects of viscosity 
contrast on buoyantly unstable miscible fluids in a porous medium, in understanding the 
effect of precipitation reactions during flow displacements in porous media in the context of 
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CO 2 sequestration techniques |T3l 133] . 
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